text = c("Memorial Day",
"Veterans Day",
"Taliban takeover Afghanistan"),
height = c(25.8,17.8,29.2))
# Generate a 5-day rolling average variable
defense_subset <- defense_subset %>%
group_by(Group) %>%
mutate(roll_mean = zoo::rollmean(Value, 5, na.pad = T))
# Rename the time series
defense_subset$Group <- recode(
as.factor(defense_subset$Group),
`national_legislators_democrat` = "Democrats in Congress",
`national_legislators_republican` = "Republicans in Congress",
`state_legislators_democrat` = "Democrat State Legislators",
`state_legislators_republican` = "Republican State Legislators")
defense_subset = subset(defense_subset, !is.na(defense_subset$roll_mean))
# Create Panel B of Figure 3
defense <- ggplot(data = defense_subset) +
geom_line(aes(x = Date, y = roll_mean*100, colour = Group,  group = Group, size = Group)) +
scale_size_manual(values = c(1.5,1.5,1,1)) +
geom_label(data = events, mapping = aes(label = text, x = time, y = height),
angle = 0,  size = 3.75) +
scale_colour_manual("Group:",
values = c(`Democrats in Congress` = "blue",
`Republicans in Congress`  = "red",
`Democrat State Legislators`  = "cornflowerblue",
`Republican State Legislators`  = "darksalmon"
)) +
ylab("Attention to defense in %") +
theme_bw() +
ylim(0,30.5)+
guides(size = "none") +
guides(colour = guide_legend(override.aes = list(size = c(3,3,2,2)))) +
theme(
axis.line = element_line(color = "black"),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 13),
axis.text.x = element_text(size = 10),
legend.title = element_text(size = 13),
legend.text = element_text(size = 9)
)
#===============================================================================
# OUTPUT
#===============================================================================
# Combine Panel A and Panel B to make Figure 3
plot_grid(immi  + theme(legend.position = "none"),
defense + theme(legend.position = "bottom",
legend.text = element_text(size = 15)),
ncol = 1, align = "v",
rel_heights = c(1,1), labels = c("A", "B"))
ggsave(paste0(root_dir, "/figures/figure01.png"), width = 14, height = 7, dpi = 300)
################################################################################
# 03-figure02.R
# Date:     Feb 6, 2025
# Authors:  Oscar Stuhler and Andreu Casas
# Purpose:  This script fits the main models for the results presented in
#           sections 5.1 and 5.2 of the paper. It then generates Figure 2.
# Note:     Models were fit on NYU's HPC cluster. Instances used 40 CPU cores
#           and 100 GB of memory. With this setup, fitting both models took
#           half a day. In this script, we comment out the portion where we fit
#           the models and instead load the fitted model objects from disk.
# Data In:
#           Two datasets in which we have group-day-issue-level information about
#           how much attention each group of analysis paid to each day to each
#           topic we study:
#           /data/group-day-issue-level-dataset-01-2018
#           /data/group-day-issue-level-dataset-01-2021
# Data Out:
#           1. Two estimated VAR models, one for 2018 and one for 2021.
#              /models/MODEL_I_2018.Rdata
#              /models/MODEL_I_2021.Rdata
#           2. Two IRF objects with the impulse response estimates based on the
#              fitted models.
#              /models/MODEL_I_IRFs_2018.Rdata
#              /models/MODEL_I_IRFs_2021.Rdata
#
#           3. Figure 2 of the main paper:
#              - ./figures/figure02.png
################################################################################
#===============================================================================
# PACKAGES
#===============================================================================
library(dplyr)
library(tidyr)
library(vars)
library(ggplot2)
library(ggh4x)
library(boot)
#===============================================================================
# Directories
#===============================================================================
local = T
root_dir = if(local){
"~/Google Drive/My Drive/Leads_follows_revision/lead_follow_240511"}else{
"/scratch/oms279/lead_follow_240511"
}
#===============================================================================
# FITTING MODELS
#===============================================================================
#
####### Defines years variable
#years = c("2021", "2018")
#
####### For each year, run the models with the same parameters
#for(year in years){
#
#  cat(year, "\n")
#  #-----------------------------------------------------------------------------
#  # DATA
#  #-----------------------------------------------------------------------------
#  # Model data
#  load_file = paste0(root_dir, "/data/group-day-issue-level-dataset-01-",
#                     year,
#                     ".csv")
#  maindb <- read.csv(load_file)
#
#  # Sort:
#  maindb = maindb[order(maindb$IssueState, maindb$Date),]
#
#  #-----------------------------------------------------------------------------
#  # MAIN
#  #-----------------------------------------------------------------------------
#  # Create a formula object that will facilitate transforming the topic variable
#  # from factor to dummies.
#  if(year == "2021"){
#    variables <- c("national_legislators", "national_media", "state_legislators",
#                   "state_media", "state_partisans", "state_random_partisans", "Trump", "IssueState")
#  }
#  if(year == "2018"){
#    variables <- c("national_legislators", "national_media", "state_legislators",
#                   "state_media", "state_partisans", "Trump", "IssueState")
#  }
#  mformula <- formula(paste0("~", paste0(variables, collapse = " + ")))
#  model_data <- model.matrix(mformula, maindb[, variables])
#  model_data <- model_data[, 2:ncol(model_data)] # remove intercept
#
#  # Splitting the covariates of interest from the issue dummy variables
#  X_endogenous <- model_data[, variables[which(!variables %in% c("IssueState"))]]
#  X_exogenous <- model_data[
#    , which(!colnames(model_data) %in%
#              variables[which(!variables %in% c("IssueState"))])]
#
#  # Define lags, number of bootstraps, and set seed
#  p <- 5
#  runs <- 300
#  seed <- 78223
#  set.seed(seed)
#
#  # Run the VAR model and estimate IRFs
#  var_model_merged <- VAR(y = X_endogenous, p = p, exogen = X_exogenous)
#  var_irfs_cum_merged <- irf(var_model_merged, n.ahead = 15, ortho = F,
#                             cumulative = TRUE, runs = runs)
#
#  # Save the estimated VAR model and IRFs.
#  save(var_model_merged, file = paste0(root_dir, "/models/MODEL_I_", year, ".Rdata"))
#  save(var_irfs_cum_merged, file = paste0(root_dir, "/models/MODEL_I_IRFs_", year, ".Rdata"))
#}
#===============================================================================
# VISUALIZE THE MODELS IN FIGURE 2 USING THE IRFs
#===============================================================================
# Load the figure code
source(paste0(root_dir,"/code-figures/00-functions.R"))
# Define which type of effect to show
which_effects <- "structural"
# Load the IRFs for Model I
year = "2021"
model_name <- paste0(root_dir, "/models/MODEL_I_IRFs_", year, ".Rdata")
print(load(model_name))
var_irfs_cum_merged_2021 = var_irfs_cum_merged
year = "2018"
model_name <- paste0(root_dir, "/models/MODEL_I_IRFs_", year, ".Rdata")
print(load(model_name))
var_irfs_cum_merged_2018 = var_irfs_cum_merged
# Generate data for plotting 5-day responses to a permanent shock
final_input_2018 <- transform_var_irf(var_irfs = var_irfs_cum_merged_2018, days = 5)
plot_db_2018 <- make_plot_data(final_input = final_input_2018,
which_effects = which_effects,
which_day = 5,
column_group = c("state_legislators"))
plot_db_2018$Year = "2018"
final_input_2021 <- transform_var_irf(var_irfs = var_irfs_cum_merged_2021, days = 5)
plot_db_2021 <- make_plot_data(final_input = final_input_2021,
which_effects = which_effects,
which_day = 5,
column_group = c("state_legislators"))
plot_db_2021$Year = "2021"
plot_db = rbind(plot_db_2018, plot_db_2021)
plot_db$rowgroup = recode(plot_db$rowgroup,  "Trump" = "President")
# Make level variable
plot_db$level <- ifelse(plot_db$rowgroup %in% c("President", "Members of\nCongress",
"National Media"),
"National\ngroup", "State\ngroup")
#  Define the order for the plot
plot_db$direction <- factor(plot_db$direction,
levels = c("Row groups' response",
"State legislators' response"
))
plot_db$year_dir = paste0(plot_db$direction, plot_db$Year, sep = "")
plot_db$rowgroup = recode(plot_db$rowgroup,
"National Media" = "National\nMedia",
"State Media" = "State\nMedia",
"State Partisans" = "State\nPartisans")
# Generate Figure
plot_db$level = factor(plot_db$level, levels = c("State\ngroup", "National\ngroup"))
plot_overall = ggplot(plot_db, aes(x = direction,  ymin = lwr,
ymax = upr, col = direction, pattern = Year)) +
geom_hline(yintercept = 0, color = "red") +
geom_linerange(aes(x=year_dir, ymin=lwr, ymax=upr), size = 5.5) +
geom_segment(aes(x=year_dir, xend=year_dir, y=pe-0.01, yend=pe+0.01), color="black", size=5) +
geom_text(aes(x=year_dir,
y = upr+.23,
label = Year), col = "black", size = 3.5) +
facet_nested(level + rowgroup ~ colgroup, switch = "y") +
coord_flip() +
ylab("\n5-day effect of a permanent 10 percentage point increase") +
scale_color_manual("Effect direction:", values = c("gray80", "gray40")) +
scale_shape_manual("Year:", values = c(2,1)) +
guides(col = guide_legend(ncol=2,nrow=1, reverse = TRUE)) +
theme_bw() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
panel.background = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank(),
strip.text.y.left = element_text(angle = 0, size = 11,
margin = margin(0,.2,0,.2, "cm")),
strip.text.x = element_text(size = 11),
axis.text.y = element_blank(),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 11),
axis.ticks.y = element_blank(),
axis.line = element_line(color = "black"),
axis.title.y = element_blank()
)
plot_overall
#  Save the figure
ggsave(paste0(root_dir, "/figures/figure02.png"), plot_overall, width = 8.5, height = 7)
################################################################################
# 06-figureD1.R
# Date:     Feb 6, 2025
# Authors:  Oscar Stuhler and Andreu Casas
# Purpose:  This script fits independent models for each issue. This provides the
#           basis for the analyses presented in Appendix G. It then generates
#           Figure G1
# Note:     Models were fit on NYU's HPC cluster. Instances used 40 CPU cores
#           and 100 GB of memory. With this setup, fitting both models took
#           half a day. In this script, we comment out the portion where we fit
#           the models and instead load the fitted model objects from disk.
# Data In:
#           Two datasets in which we have group-day-issue-level information about
#           how much attention each group of analysis paid to each day to each
#           topic we study. MoCs were differentiated by the state they, as in
#           the analyses for Figures 3 and 4.
#           are associated with.
#           /data/group-day-issue-level-dataset-04-2018.csv
#           /data/group-day-issue-level-dataset-04-2021.csv
# Data Out:
#           1. Two list objects, containing the estimated models and IRFs
#              for each issue.
#              /models/MODEL_by_issue-2018.Rdata
#              /models/MODEL_by_issue-2021.Rdata
#
#           2. Figure G1 of the Appendix:
#              - ./figures/figureG1.png
################################################################################
#===============================================================================
# PACKAGES
#===============================================================================
library(dplyr)
library(tidyr)
library(vars)
library(boot)
library(stringr)
library(ggplot2)
library(ggh4x)
library(cowplot)
#===============================================================================
# Directories
#===============================================================================
local = T
root_dir = if(local){
"~/Google Drive/My Drive/Leads_follows_revision/lead_follow_240511"}else{
"/scratch/oms279/lead_follow_240511"
}
#===============================================================================
# FITTING MODELS
#===============================================================================
#for(year in c("2018", "2021")){
#  cat(year, "\n")
#  load_file = paste0(root_dir, "/data/group-day-issue-level-dataset-04-",
#                     year,
#                     ".csv")
#
#  maindb <- read.csv(load_file)
#
#  # Sort order:
#  maindb = maindb[order(maindb$IssueState, maindb$Date),]
#
#  #-------------------------------------------------------------------------------
#  # MAIN
#  #-------------------------------------------------------------------------------
#  # Define lags, number of bootstraps, and set seed
#  p <- 5
#  runs <- 300
#  seed <- 78223
#  set.seed(seed)
#
#  # Generate an issue domain variable
#  if(year == "2021"){
#    variables <- c("national_legislators", "national_media", "state_legislators",
#                   "state_media", "state_partisans", "state_random_partisans", "Trump", "IssueState")
#  }
#  if(year == "2018"){
#    variables <- c("national_legislators", "national_media", "state_legislators",
#                   "state_media", "state_partisans", "Trump", "IssueState")
#  }
#
#  # Run three models for the different types of issues
#  issue_results <- list()
#  unique_issues = unique(str_remove(maindb$IssueState, "[A-Z]{2}$"))
#
#  for(issue_i in unique_issues){
#    cat(issue_i, "\n")
#
#    # Subset the data by state group
#    issue_db <- maindb[which(str_remove(maindb$IssueState, "[A-Z]{2}$") == issue_i),]
#
#    mean_vec = c()
#    for(group_i in variables[1:7]){
#      issue_results[[issue_i]][[group_i]] = exp(mean(issue_db[[group_i]]))
#    }
#
#
#    # Generate formula object
#    mformula <- formula(paste0("~", paste0(variables, collapse = " + ")))
#    model_data <- model.matrix(mformula, issue_db[, variables])
#    model_data <- model_data[, 2:ncol(model_data)] # removing intercept
#
#    # Splitting the covariates of interest from the issue dummy variables
#    X_endogenous <- model_data[, variables[which(!variables %in% c("IssueState"))]]
#    X_exogenous <- model_data[
#      , which(!colnames(model_data) %in%
#                variables[which(!variables %in% c("IssueState"))])]
#
#    # Run the VAR model and estimate IRFs
#    var_model_merged <- VAR(y = X_endogenous, p = p, exogen = X_exogenous)
#    var_irfs_cum_merged <- irf(var_model_merged, n.ahead = 15, ortho = F,
#                               cumulative = TRUE, runs = runs)
#
#    # Store results in list in list
#    issue_results[[issue_i]][["IRF_res"]] = var_irfs_cum_merged
#  }
#
#  # Save the estimated VAR model and IRFs.
#  object_name = paste0(root_dir, "/models/MODEL_by_issue-", year, ".Rdata")
#  save(issue_results, file = object_name)
#}
#===============================================================================
# VISUALIZE THE MODELS IN FIGURE G1 USING THE IRFs
#===============================================================================
# Load utils & define which effects to show
source(paste0(root_dir,"/code-figures/00-functions.R"))
which_effects <- "structural"
plot_db = data.frame()
for(year in c("2018", "2021")){
##### Load
cat(year, "\n")
object_name = paste0(root_dir,
"/models/MODEL_by_issue-", year, ".Rdata")
print(load(object_name))
for(issue_i in names(issue_results)){
cat(issue_i, "\n")
results <- transform_var_irf(var_irfs = issue_results[[issue_i]]$IRF_res, days = 5)
plot_db_i <- make_plot_data(final_input = results,
which_day = 5,
which_effects = which_effects, column_group = c("state_legislators"))
plot_db_i$issue = issue_i
plot_db_i$Year = year
for(group_i in names(issue_results[[issue_i]])[1:5]){
plot_db_i[[group_i]] = issue_results[[issue_i]][[group_i]]
}
plot_db = rbind(plot_db, plot_db_i)
}
}
# Generate Figure 6
p1 = ggplot(subset(plot_db, plot_db$rowgroup == "Members of\nCongress" & plot_db$direction == "State legislators' response")) +
geom_point(aes(x = national_legislators*100, y = pe,
col = Year)) +
facet_nested(.~"State legislators' response to Members of Congress") +
stat_smooth(aes(x = national_legislators*100, y = pe, group = Year, col = Year),
method = "lm", se = F
#fill = "grey90", alpha = .5
) +
xlab("Issue prevalence in percent among\nMembers of Congress") +
ylab("5-day effect of a permanent 10\npercentage point increase") +
scale_color_manual("Year:",
breaks = c("2021", "2018"),
values = c("gray80", "gray40")) +
#scale_shape_manual("Year:", values = c(2,1)) +
guides(col = guide_legend(reverse = T)) +
theme_bw() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 11),
axis.ticks.y = element_blank(),
axis.line = element_line(color = "black")
)
# Generate Figure G1
p2 = ggplot(subset(plot_db, plot_db$rowgroup == "Members of\nCongress" & plot_db$direction == "Row groups' response")) +
geom_point(aes(x = national_legislators*100, y = pe,
col = Year)) +
facet_nested(.~"Members of Congress' response to State Legislators") +
stat_smooth(aes(x = national_legislators*100, y = pe, group = Year, col = Year),
method = "lm", se = F
#fill = "grey90", alpha = .5
) +
xlab("Issue prevalence in percent among\nMembers of Congress") +
ylab("5-day effect of a permanent 10\npercentage point increase") +
scale_color_manual("Year:",
breaks = c("2021", "2018"),
values = c("gray80", "gray40")) +
guides(col = guide_legend(reverse = T)) +
theme_bw() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 11),
axis.ticks.y = element_blank(),
axis.line = element_line(color = "black")
)
p3 = ggplot(subset(plot_db, plot_db$rowgroup == "Members of\nCongress" & plot_db$direction == "State legislators' response")) +
geom_point(aes(x = state_legislators*100, y = pe,
col = Year)) +
facet_nested(.~"State legislators' response to Members of Congress") +
stat_smooth(aes(x = state_legislators*100, y = pe, group = Year, col = Year),
method = "lm", se = F
#fill = "grey90", alpha = .5
) +
xlab("Issue prevalence in percent among\nState Legislators") +
ylab("5-day effect of a permanent 10\npercentage point increase") +
scale_color_manual("Year:",
breaks = c("2021", "2018"),
values = c("gray80", "gray40")) +
guides(col = guide_legend(reverse = T)) +
theme_bw() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 11),
axis.ticks.y = element_blank(),
axis.line = element_line(color = "black")
)
# Generate Figure G1
p4 = ggplot(subset(plot_db, plot_db$rowgroup == "Members of\nCongress" & plot_db$direction == "Row groups' response")) +
geom_point(aes(x = state_legislators*100, y = pe,
col = Year)) +
facet_nested(.~"Members of Congress' response to State Legislators") +
stat_smooth(aes(x = state_legislators*100, y = pe, group = Year, col = Year),
method = "lm", se = F
) +
xlab("Issue prevalence in percent among\nState Legislators") +
ylab("5-day effect of a permanent 10\npercentage point increase") +
scale_color_manual("Year:",
breaks = c("2021", "2018"),
values = c("gray80", "gray40")) +
guides(col = guide_legend(reverse = T)) +
theme_bw() +
theme(
legend.position = "bottom",
legend.text = element_text(size = 12),
legend.title = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 11),
axis.ticks.y = element_blank(),
axis.line = element_line(color = "black")
)
by_plots = plot_grid(p1 + theme(legend.position = "none",
plot.margin = margin(.75,.38,.1,.1, "cm"),
axis.title.x = element_blank(),
axis.text.x = element_blank()),
p3 + theme(legend.position = "none",
plot.margin = margin(.75,.38,.1,.1, "cm"),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank()),
p2 + theme(legend.position = "none",
plot.margin = margin(.75,.38,.1,.1, "cm")),
p4 + theme(legend.position = "none",
plot.margin = margin(.75,.38,.1,.1, "cm"),
axis.title.y = element_blank(),
axis.text.y = element_blank()),
rel_heights = c(1,1.2), rel_widths = c(1.13,1),
labels = c("A", "B", "C", "D"),
ncol = 2, nrow = 2)
by_issue_p = plot_grid(by_plots,
get_legend(p4),
nrow = 2,
rel_heights = c(1,.1))
# Save the Figure G1
plot_name = paste0(root_dir, "/figures/figureG1", ".png")
ggsave(plot_name, by_issue_p, width = 7.5, height = 7.5)
